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The destruction of anomalous diffusion of the Harper model at criticality, due to weak nonlin- 
earity x> is analyzed. It is shown that the second moment grows subdiffusively as {m.2) ~ t" up 
to time t* ~ x 1 ■ The exponents a and 7 reflect the multifractal properties of the spectra and the 
eigenfunctions of the linear model. For t > t* , the anomalous diffusion law is recovered, although 
the evolving profile has a different shape than in the linear case. These results are applicable in 
wave propagation through nonlinear waveguide arrays and transport of Bose-Einstein condensates 
in optical lattices. 
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I. INTRODUCTION 

The spreading of a quantum mechanical wavepacket 
for a particle moving in a periodic lattice is a textbook 
example. One finds that the width increases ballistically 
with time and the corresponding cigenstatcs are extended 
(Bloch states). In the opposite case of strongly disor- 
dered systems, the eigenstates are exponentially local- 
ized, resulting in a total halt of the wavepacket spreading 
in the long-time limits. Between these extremes flour- 
ishes the world of quantum systems with anomalous dif- 
fusion. These include systems studied in the early days 
of quantum mechanics, such as Bloch electrons in a mag- 
netic field^ as well as quasicrystals^ and disordered sys- 
tems at the metal-insulator transition^. In these cases, 
the eigenfunctions (or even the spectrum) show a fractal 
structure which dictates the wavepacket spreading. 

The physical motivation to study such systems, along 
with the mathematically intriguing nature of their spec- 
tra, led to a series of works that eventually advanced our 
understanding of their dynamical properties. Most of 
these works have focused on the analysis of the temporal 
decay of the survival probability P(t) at the initial posi- 
tion no and the growth of the wavepacket second moment 

Despite all this activity, nothing is known about the ef- 
fect of nonlincarity in the wavepacket spreading for sys- 
tems with fractal spectra and eigenfunctions. Notice- 
able exceptions are Refs. 9 and 10 which numerically 
studied the wavepacket dynamics for the prototype one- 
dimensional (ID) tight-binding Harper model with non- 
linearity. Their conclusions, however, are contradictory. 
Although both studies conclude that infinitesimal non- 
lincarity results in a short-time subdiffusive spreading 
of the variance 777,2 ~ £ Q , they give different values for 
the power-law exponent. While Rcf. 10 concludes that 
the subdiffusive spreading persists for longer time, Ref. 9 
reports a saturation of the second moment. Moreover, 
they both fail to identify traces of the underlying linear 
model's mult ifr act ality in the dynamics generated by the 
corresponding nonlinear one. 

In the present paper, we address the nonlincarity- 



induced destruction of the anomalous diffusion. Our fo- 
cus is on the ID Harper model at criticality but we expect 
that our results will be applicable for other critical mod- 
els (such as Fibonacci lattices) as well. Our detailed nu- 
merical calculations and theoretical arguments show that 
there is a critical nonlincarity x* , above which the nature 
of anomalous diffusion (associated with the linear model) 
is altered. In particular, we found that for paramctrically 
large time, the average (over initial sites) spreading of an 
initially localized 5— like wavepacket spreads as 



(m 2 (t))(xt a with 2(3 



£>2 < a < 2/3 



(1) 



where D% is the the fractal dimension of the Local Den- 
sity of States (LDoS) and a is independent of the non- 
lincarity x- The time scale t* up to which Eq. fTJ is 
observed depends on x as 



t* 



X 1 



with 7 = 



(2) 



Beyond this time scale, the effects of nonlinearity dimin- 
ish and we recover the spreading of the linear model i.e. 
(ni2(t)) cx i 2/3 . However, other moments still behave dif- 
ferently with respect to the linear (i.e. X = 0) case. Fi- 
nally, we show that for any x 7^ 0j the core of the evolving 
profile satisfies the scaling behavior, 



P s (x,t) = y/(m 2 (t))P x (n, f ) , x 



n 



(3) 



This paper is structured in the following way. In Sec. 

II an overview of the linear Harper model is presented. 
The nonlinear Harper model will be introduced in Sec. 

III and the effects of nonlinearity on the spreading and 
the evolving profile will be investigated both numerically 
and analytically. Finally, wc will draw our conclusions in 
Sec. IV. 



II. THE HARPER MODEL: OVERVIEW 

In this section, we will briefly review the basic proper- 
ties of the linear Harper Model. Its dynamics is described 



2 



by the standard ID tight-binding model, 

MM*) 



dt 



ip n+1 (t) + ip n ^(t) + V n ip n (t) 



(4) 



Above, ip n (t) denotes the probability amplitude for a par- 
ticle to be at site n at time t and we will assume that the 
on-site potential V n takes the form: 



V n = A C0S(27T<771 



(5) 



When a is an irrational number, the period of the on-site 
potential V n is incommensurate with the lattice period. 
In such cases, the eigenstates of the linear system are 
extended when A < 2 and the spectrum consists of bands 
(ballistic regime). For A > 2, the spectrum is point- 
like and all states are exponentially localized (localized 
regime). The most interesting case is the critical point 
A = 2, where we have a metal-insulator transition. At 
this point, the spectrum is a zero measure Cantor sefci 1 -, 
while the eigenstates show self-similar fluctuations on all 
scale a 12 ! 13 . 

For the critical case, where A = 2, previous studies 
show that the temporal decay of the probability to remain 
at the original site n up to time t goes as 



P(t)^\^ no (t)\ 2 ~l/t D * 
while the wavepacket second moment grows as 



(G) 



m 2 (t) = Y,(n-nof\Mt)\ 2 ~ t 20 with [3 > D% / Df 

n 

(7) 

where the exponent f3 depends on both D% and the 
correlation dimension of the fractal eigenfunctions £ 
In fact, recent studies were able to demonstrate that 
the center of the wavepacket spatially decreases as \n — 

nol ^ 1 for | n — no\ <C t 13 while the front shape scales 
as£ 

P(n,t) = A(t) exp(-|(n - n o )/V^ 2 \ 1/(1 - 0) ) (8) 
where A(t) is the height of the wave profile at time t. 



III. THE NONLINEAR HARPER MODEL 

Introducing an on-site nonlinear term into Eq. ([5]), 
we come out with the nonlinear Harper model (NLHM). 
Mathematically, it is described by the following discrete 
nonlinear Schrodinger (DNLS) equation, 



1 dt 



VWl (*) + (*) + K (*) ~ Xl VVi (*) I Vn (*) 

(9) 

We want to investigate the temporal behavior of the 
variance (m2(i))0 for the NLHM at the critical point 
A = 2. The initial condition is always taken to be a 
(5-likc localized state i.e. tp n (t = 0) = S n no . An aver- 
aging (-)^i over different initial positions hq (at least 50) 



of the S function [or cquivalently over a random distri- 
bution of the phase (j> in Eq. ([5])] is done. In our calcu- 
lations, we assume a to be equal to the golden mean 
o"G = (v5 — l)/2. Equation §§§ has been integrated 
numerically using a finite-time-step fourth order Rungc- 
Kutta algorithm on a self-expanding lattice in order to 
eliminate finite-size effects^ 4 -. Whenever the probability 
of finding the particle at the edges of the chain exceeded 
ICT 10 , ten new sites were added to each edge. Numer- 
ical precision is checked by monitoring the conservation 
of probability (norm) ^2 n \ip n (t)\ 2 = 1. In all cases, the 
deviation from unity was less than 10 -6 . 



A. Wavepacket Spreading 

In Fig. [1] we report our numerical results for some 
representative \ values in a double logarithmic plot. In 
all cases, the variance displays a power-law behavior 
(m2(t)) ~ t a . Leaving aside the very initial spreading 
t < 1 (which in all cases is ballistic), we observe that 
for x > x* ~ 5.5 x 10 -4 the spreading is subdiffusive 
with a < 2(3. This spreading is valid up to time t < t* 
and then relaxes to the anomalous diffusion a — 2/3 that 
characterizes the linear Harper model. 

To obtain the law of spreading of Eq. (Q]) for t < t*, 
we consider the following simple model: The initial site 
together with the nearby sites n = n ±5n is considered as 
a confined source 1 ^. Anything emitted from them moves 
according to the spreading law of the linear Harper i.e. 
s(t) = (rri2{t)) < t > = vt@ since for t < t* the leaking 
norm is small and therefore the nonlinear term in Eq. (J9]) 
is negligible with respect to the on-site potential V n . 

Initially, all probability is concentrated at the source. 
Due to the nonlinear nature of the system, the decay 
rate is not constant but depends on the remaining norm. 
Thus, the decay process is characterized by the nonlinear 
equation, 



dP{t)/dt = -T P(t) P(t) 



(10) 



which leads to non-exponential decay^ On the other 
hand, for x — 0- the decay is power-law 



P{t) - l/t s 



(11) 



with S = Z?2 [sec Eq. ^ above]. It is, therefore, natural 
to expect that for x 7^ 0, we will have a slower decay 
with S < £>2 due to self-trapping] 17 i 18 i 19 The variance of 
the confined-source model is then given by 2 ^ 

/■OO r-t 

M PS (i)« / dss 2 dt'(-P(t'))6(s-v(t-t'f), (12) 
Jo Jo 

where —P(t) is the flux emitted from the confined source. 
Substituting Eq. {TED for P(t) we get 



Mps(t) - v 1 



dt( 7 ) 
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(t-t) 



2D 



(13) 
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FIG. 1: (color online) The variance (7712(^)0 for various x 
values of the NLHM. The \ = 10 -4 case is shifted downwards 
in order to distinguished it from the x = case. Dashed lines 
have slopes as indicated in the figure and are drawn to guide 
the eye. Inset: the decay of P(t) ~ l/t^ 2 for the linear 
Harper model (x = 0). 



FIG. 2: The fitting values (stars) of the power-law exponent 
a vs. V for the Fibonacci model. Circles are the fitting expo- 
nent 2/3 of the linear model for the spreading (matt)) ~ t 2 ' 3 , 
while diamonds are the extracted exponents 2/3 — D%, where 
D2 was obtained from the decay of P(t) ~ l/t D 2 of the linear 
model. 



If we do a variable substitution t = t /t, we will get 
Mps(i) ~ v 2 t 20 - s f dt (l/f) m (l - f) 2/3 , (14) 

JQ 

which leads to Eq. (JU with a = 2/3 - 6 > 2/3 - D%. 

In order to compare the results of the numerical simu- 
lations with the theoretical predictions of Eq. (p} we have 
calculated for the linear Harper model the decay expo- 
nent (D^)^ of the survival probability (see inset of Fig. |T|) 
and the power-law exponent of the variance (777,2 (£))</,. 
The least-squares fit gives the values (D^)^ ~ 0.30 ±0.03 
and 2/3 « 1.00±0.03. The numerically extracted value of 
a = 0.71 fulfills nicely the bound 2(3 - D% = 0.70 given 
by Eq. ©. 

To further test the validity of Eq. JT]), we have also 
performed simulations with the Fibonacci model, where 
£>2 an d D% (and, therefore, /3) can be varied according to 
the on-site potential V n . The latter takes only two values 
±V (V 0) that are arranged in a Fibonacci sequence.^ 
Again, we find a power-law spreading (7712(f)) ~ t a . 
The extracted exponents a corresponding to various V's 
are reported together with the theoretical predictions in 
Fig. [2] and they confirm the validity of Eq. (fTJ) . 

From Fig. [T] we see that the time scale t* up to 
which the power law of Eq. (fTJ) applies depends on the 
strength of the nonlincarity parameter \- One can esti- 
mate the time t* from the fact that the effective potential 
V x = — xlV'noCi)! 2 is comparable with the on-site poten- 
tial V no ~ A of the linear model at t*. After this time, one 
expects that the effect of nonlinearity on the wavepacket 
spreading is negligible and, therefore, the survival prob- 
ability should decay as P(t) = |^„ (f)| 2 ~ l/t D *. Fol- 
lowing this line of argumentation, we have that 

v no - xlVUOl 2 - v no ~ x0-/t*) D Z, (15) 




t/t*(x) 



FIG. 3: (color online) The variance vs. scaled time, t* is 
extracted by scaling the time-axis so that the variance curves 
overlap in the time-region where Eq.JT} is valid. Inset: the 
scaling of t* versus x- The dashed line indicates the least 
square fit with slope 3.96 [the theoretical prediction Eq. ([2]) 
estimates it to be 3.45 ± 0.35]. 



leading to Eq. ([2]). To test this theoretical prediction, wc 
have manually scaled the time axis so that the variance 
curve where Eq. |T]) applies overlaps for various x values. 
The extracted scaling parameters t* (x) are plotted in the 
inset of Fig. [3] One can clearly see that the numerical 
data confirms the theoretical prediction of Eq. (J5J). 

Finally, we will discuss the shape of the evolving wave- 
function for x > X* ■ Although for t > t* the temporal 
behavior of the variance becomes the same as that of the 
linear Harper model, other moments differ. The profiles 
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FIG. 4: (color online) The scaled probability distributions 
P 3 (x) versus x of the NLHM at different time t (including 
t > t*) and for various \ > X* ■ The nice overlap of the curves 
at the core |a;| < 8 confirms the validity of the scaling law 
Eq. 0. The scaling is lost for \x\ > 8. The lower inset shows 
the behavior near the origin. The blue solid line is the fitting 
curve Eq. Q5) . For comparison we also report in the upper 
inset the probability distribution P s (x) of the linear Harper 
model. 

reported in Fig. U are snapshots of the average wavefunc- 
tion at various times and for three representative values 
of x = 1,3 and 5 (the linear case \ = is also plotted in 
the upper inset for comparison). They are plotted with 
the scaling assumption in Eq. The data collapse for 
| a; | < 8 (different curves corresponding to various times 
t and nonlinearity strengths x) reveals that this repre- 
sentation is not affected much by finite- x corrections, al- 
though the tails of the profile are x-dependent. In fact, 
we found that the probability distribution near the center 
can be described by the formula: 

P s (aO ~M~ 7 exp(-|x|Ax>), M<8, (16) 

where 7 w 0.7 and 1.82. We note that a similar 

expression for the core of the probability distribution ap- 
plies for ID and quasi-lD disordered models with zero 
nonlinearityi 14 ' 27 



in Fig. [5] We see that up to x* ~ 5.5 x 10~ 4 , the time 
average survival probability (P(T))x remains unchanged. 
As the nonlinearity strength x is increased further, a frac- 
tion of the excitation begins to localize at the initial site. 
As a result, the fraction of the excitation that can prop- 
agate is now effectively smaller, leading to smaller values 
of (rri2(t)) as x is increased. We note that a similar type 
of self-trapping phenomeno n 17 ! 18 was observed in various 
nonlinear lattices j£ i 10 ' 19 i 21 ' 22 i 23 ' 24 i 25 though in all these 
cases the value of x* was much larger, i.e. x* ~ 0(1). 

The following heuristic argument provides some under- 
standing of the appearance of self-trapping phenomenon 
for the NLHM. We consider successive rational approxi- 
mants cr.; = Pi/qi of the continued fraction expansion of 
a. On a length scale the periodicity of the potential is 
not manifested and we may assume that for x < X* t ne 
cigenfunctions preserve their critical structure as in the 
case of x — 0. In this case, the partitioning of the energy 
over the qi sites is^ 

n = -f E iv^+E^n-i+C-tVvo+E Vn i^™i 2 

71 — 1 71— 1 71 — 1 

(18) 

Let us first estimate the energy Tt ss associated with the 
initial state 5 n ^ no . Since the probability distribution is 
located at a single site, we get an energy 

(W->* = -f (19) 

Now we want to evaluate the partitioning of the en- 
ergy 7i ox t for a fractal wavefunction. In this case, the 
first term in Eq. (|18p is the inverse participation num- 
ber, which scales with the system size qi as 

P2 ~ q ; Dt (20) 

Hence, in the thermodynamic limit where qi — > 00, this 
term will go to zero. We define the sum of the two other 
terms as follows: 

n 

S = <]T(CVn-l + C-lV-n + V n |^„| 2 )) (21) 
n—l 



B. Critical Nonlinearity 

Going back to Fig. [TJ we observe that the destruction 
of the anomalous diffusion of the linear model takes place 
for x > X* ■ To quantify x*> we evaluate the time- average 
survival probability (P(T))t, defined asi 8 .^^ 

(P(T)) T = lim i f \^ no {t)\ 2 dt (17) 

for various values of the nonlinearity strength x- In our 
numerical calculations, we took an average over a time 
interval of T = 20000. Our numerical results are reported 



We found from numerical calculations (sec inset of Fig. [5]) 
that S goes to a finite value in the thermodynamic limit, 
and, therefore, arrive at the equation {TL e ^t)<t> = 

If Ticxt is higher than the initial energy H ss , then con- 
servation of energy prevents the partitioning of energy 
over all qi sites. Therefore, in the thermodymanic limit, 
we get for the critical nonlinearity x* that 

X* ~ -25. (22) 

The numerical data in Fig. [5] (inset) indicate that S ~ 
—2.2 x 10 -4 , which is consistent with our numerical eval- 
uation of x* ~ 5-5 x 10~ 4 . 
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FIG. 5: The integrated survival probability (P(T))t vs. \. 
For x > 5.5 x 10~ 4 we observe an increase in (P(T))t- The 
first point in the curve corresponds to the case x — an d is 
included in the log-log plot as a reference point. Upper inset 
shows the scaling of the P2 and —S vs. the system size qi. 
The horizontal dashed line is S = —2.2 x 10~ 4 . 



nent of the temporal spreading of the wavcpacket vari- 
ance. These bounds reflect the fractal dimension of the 
LDoS of the linear system. This nonlinear spreading ap- 
pears for nonlinearity strength above some value and per- 
sists up to time t* ~ x 1 ^ 2 1 which depends parametri- 
cally on the nonlinearity strength. After this time, the 
linear spreading of the second moment is restored; other 
moments, however, are still affected by the nonlinearity. 
For the central part of the evolving profile, we have also 
found a scaling relation that applies to any time and any 
nonlinearity strength. 



The importance of these results lies in their applicabil- 
ity to quasipcriodic photonic structures (such as optical 
super-lattices^) and arrays of magnetic micro-traps for 
atomic Bose-Einstein Condensates^ (BEC). In the lat- 
ter case, nonlinear wave self-action appears due to atom- 
atom scattering in BEC, while in optical crystals, it ap- 
pears due to light-matter interactions. 



IV. CONCLUSIONS 

In this paper, we studied the nonlinear Harper model 
at criticality and found bounds for the power-law expo- 
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